REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMQNo .  0704-0188 

The  public  reporting  burden  for  this  collection  of  information  IS  estimated  \Q  average  1  hour  par  response,  including  the  Tims  fof  reviewing  in  si  rue  Irani,  searching  emsTing  data  sources, 
gal  he  ring  and  maintaining  th*  data  Jifrtdad,  and  completing  and  reviewing  The  collection  of  information.  Send  comments  regarding  this  burden  estimate  Of  any  Other  aspect  nl  this  collection  of 
information.  including  suggestions  for  reducing  the  burden,  to  the  Department  of  Defense,  Executive  Services  and  Communications  Directorate  (0704  0180)  Respondents  should  be  aware 
that  notwithstanding  any  other  provision  of  law.  no  person  shall  be  subject  to  any  panalty  for  failing  to  comply  with  a  collection  of  information  if  H  dots  not  display  a  currently  valid  OMB 
control  number 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ORGANIZATION. 

1.  REPORT  DATE  (OD-MM-YYYY)  2  REPORT  TYRE 

16-03-2007  Final  Report 

3.  DATES  COVERED  (From  ■  To) 

1 5  Jul  2006  -  30  Nov  2006 

4,  TITLE  AND  SUBTITLE 

Computational  Investigation  of  Unsteady  Low-Reynolds  Number 
Aerodynamics  for  Micro  Air  Vehicles 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

FA9550-06- 1-0491 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

Wei  Shyy 

Yongsheng  Lian 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME[S]  AND  ADDRESS(ES1 

Department  of  Aerospace  Engineering 

University  of  Michigan,  Ann  Arbor,  Ml  48109 

B.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSOR  INGrMONITORING  AGENCY  NAMEISI  AND  ADDRESSIES1 

AFOSR/NA 

875  N.  Randolph  Street,  Suite  325,  Room  3112 

Arlington,  VA  22203- 176jS 

Qhc+4r 

10.  SPONSORIMONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 

NUMBER(S) 

12.  DISTRfBUTlON/A  VAIL  ABILITY  STATEMENT 

Distribution  Statement  A:  Distribution  Unlimited;  Public  Release  AFRL-SR- AR-TR-07-G220 


13.  SUPPLEMENTARY  NOTES 


14,  ABSTRACT 

A  Navier-Stokcs  equation  solver  is  utilized  to  investigate  the  aerodynamics  of  a  flapping  airfoil.  The  roles  of  plunging  and  pitching 
amplitude  and  frequency,  and  Strouhal  number  are  studied.  For  a  NACA001 2  airfoil  at  zero  geometric  angle  of  attack,  chord 
Reynolds  number  of  20,000  and  at  a  given  plunging  frequency,  cither  drag  or  thrust  is  produced,  depending  on  the  plunging 
amplitude.  When  drag  is  produced,  the  viscous  force  dominates  the  total  drag  with  decreasing  influence  as  the  plunging  amplitude 
increases.  At  the  considered  plunging  amplitude  (from  0.0125c  to  0.075c),  flow  history  has  more  influence  than  the  kinematic  angle 
of  attack  in  determining  lift.  For  an  airfoil  experiencing  combined  plunge  and  pitch  motion,  both  thrust  and  input  power  increase 
with  the  Strouhal  number  (within  the  range  of  0.03  to  0.5).  For  the  case  studied,  the  thrust  is  induced  by  the  lift,  which  follows 
approximately  the  curve  of  kinematic  angle  of  attack.  Leading  edge  vortex  moves  downstream  and  interacts  with  the  trailing  edge 
vortex.  The  impact  of  a  gust  on  a  stationary  and  flapping  airfoil  is  also  studied. 

15.  SUBJECT  TERMS 

Flapping  airfoil,  pitch  and  plunge,  Low  Reynolds  number,  NACA  001 2  airfoil,  CFD 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

ACCTDATT 

18.  NUMBER 

nc 

19a.  NAME  OF  RESPONSIBLE  PERSON 

o.  REPORT 

b.  ABSTRACT 

c.  THIS  PAGE 

Abb  1  t 

Ur 

PAGES 

20 

U 

U 

u 

u 

19b.  TELEPHONE  NUMBER  (Include  area  code J 

Standard  Form  298  [Rev.  8/98) 

Prescribed  by  ANSI  Sid  239  10 
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Final  Report:  FA9550-06-1-049I 

Wei  Shyy  and  Yongsheng  Lian 
Department  of  Aerospace  Engineering 
University  of  Michigan,  Ann  Arbor,  MI  48109 

Project  Summary 

Jl  is  known  that  plunging  airfoil  can  produce  both  lift  and  thrust  with  certain  combination  of  plunging 
amplitude  and  frequency.  Motivated  by  our  interest  in  micro  air  vehicles  (MAVs),  we  utilize  a  Navier- 
Stokes  equation  solver  to  investigate  the  aerodynamics  of  a  flapping  airfoil.  The  roles  of  the  plunging  and 
pitching  amplitude  and  frequency,  and  Strouhal  number  are  studied.  For  a  symmetric  plunging  airfoil 
NACA0012  at  zero  geometric  angle  of  attack  and  chord  Reynolds  number  of  2*104s  at  the  same  plunging 
frequency,  it  can  produce  either  drag  or  thrust  depending  on  the  plunging  amplitude.  At  the  considered 
plunging  amplitude  (from  0.0125c  to  0.075c),  the  How  history  has  more  influence  than  the  kinematic  angle 
of  attack  to  determine  the  lift.  When  drag  is  produced,  the  viscous  force  dominates  the  total  drag  with 
decreasing  influence  as  the  plunging  amplitude  increases.  For  an  airfoil  experiencing  combined  plunge  and 
pitch  motion,  both  thrust  and  input  power  increase  with  the  Strouhal  number  (within  the  range  of  0,03  to 
0.5).  For  the  case  studied,  the  thrust  is  induced  by  the  lift,  which  approximately  follows  the  curve  of  the 
kinematic  angle  of  attack.  Leading  edge  vortex  moves  downstream  and  interacts  with  the  trailing  edge 
vortex.  We  also  study  the  impact  of  gust  on  stationary  airfoil  and  flapping  airfoil.  Within  the  range  of  the 
parameters  tested,  for  stationary  airfoil  the  lift  is  in  phase  with  the  velocity  but  the  drag  is  slightly  out  of 
phase.  For  flapping  airfoil,  neither  lift  nor  drag  is  in  phase  with  the  velocity. 

Nomenclature 

C})  =Drag  coefficient  per  unit  span 

Q  =Lift  coefficient  per  unit  span 

C/»  -input  power  coefficient 

Cfjnun  time-averaged  input  power  coefficient 

Cf  =thntst  coefficient 

Cfj Mm„  ^ime-averaged  thrust  coefficient 

c  =Chord  length 

/  -frequency 

hn  =non-dimenstonal  plunge  amplitude 

k  ^reduced  frequency,  &)c/2U 

L  -lift  per  unit  span 

M  ^-pitching  moment 

P  =input  power  per  unit  span 

St  -Strouhal  number,  2 fch^  IV  -  2 kh^  /  n 
T  -thrust  per  unit  span,  time  period  per  cycle 
t  =time 

Uti  ^reference  velocity,  ty  pically  freestream  velocity 
a  ^kinematic  angle  of  attack 

a{)  -  nom  i  n  a  I  ang  I  e  of  attac  k 

=effective  angle  of  attack 
-geometrical  angle  of  attack 

/?  -kinematic  angle  of  attack  due  to  pitching  motion 

tj  ^propulsive  efficiency 

&  =p  itching  angle  of  attack 

tp  =phase  difference  between  pitching  and  plunging 

fo  =circular  frequency 

Ini  induction 


Micro  Air  Vehicles  (MAVs),  which  typically  has  a  dimension  of  less  than  15  cm  and  flight  speed  around 
10  m/s,  operate  at  Reynolds  numbers,  based  on  the  mean  airfoil  chord,  of  1 50,000  and  lower.  In  some  parts 
of  the  flight  envelop  the  Reynolds  number  can  be  lower  than  I04.  The  success  of  MAV  design  requires 
collaborations  across  different  disciplines  including  aerospace  engineering,  biology,  and  electronic 
engineering.  From  an  aerodynamic  point  of  view,  two  prominent  issues  needed  to  be  understood.  First, 
how  to  generate  sufficient  lift  to  keep  the  vehicles  airborne;  second,  how  to  maintain  a  stable  flight  under 
gusty  flight  conditions.  Due  to  their  low  inertia  (2.7  g  to  O(l02)  g)  and  slow  flight  speed  (0—20  m/s)  the 
MA  Vs  are  sensitive  to  the  wind  gust. 

For  lift  generation,  researchers  have  tried  three  ty  pes  of  wings,  namely,  fixed  wing,  rotary  wing,  and 
flapping  wing.  Fixed  wings  have  been  demonstrated  in  MAVs  with  a  maximal  dimension  around  10-15 
cm.  Figure  1(a)  shows  the  15-cm  MAV  design  by  Ifju  ct  al  at  the  University  of  Florida,  However,  a 
fixed  wing  MAV  does  not  possess  the  ability  to  hover  in  a  controlled,  efficient  manner,  which  is  important 
to  many  MAV  missions.  The  rotary  wing  MAV,  namely,  a  miniature  helicopter,  can  perform  controlled 
hovering.  To  date,  a  vehicle  with  6  cm  rotor  has  been  successfully  designed*3*  Figure  1(b)  shows  a  rotary 
wing  MAV  with  6  cm  rotary  diameter.  Flapping  flight  vehicle,  aimed  at  simultaneously  integrating  lift, 
propulsion  and  control*  offers  substantial  potential  for  forward  flight,  hovering,  and  wind  gust 
adaptation.*4,5'*1 

There  is  a  substantial  gap  between  the  vehicle  performance  goals  and  the  design  capability.  In  this  paper, 
we  investigate  the  plunging  airfoil  aerodynamics,  especially  in  regard  to  thrust  and  drag  characteristics.  In 
thrust  generation,  the  published  studies  have  focused  on  the  Reynolds  number  in  the  range  between  6*  1G1 
to  4* | o4.0Kyi,,l  It  is  noted  that  the  flapping  wing  MAVs  can  employ  either  a  fully  integrated  design 
utilizing  flapping  to  generate  lift  as  well  as  thrust,  or  a  hybrid  design  utilizing  flapping  for  thrust  while 
relying  on  fixed  wings  for  lift.  Figure  1(c)  is  a  biplane  MAV  design  by  Jones  and  Platzer**1,  with  a  25  cm 
span.  In  this  design,  the  flapping  wings  generate  thrust  while  the  stationary  wing  provides  lift.  The  biplane 
design  makes  the  vehicle  less  sensitive  to  gust  and  more  resistant  to  stall.**1  The  present  paper  is  closely 
related  to  flapping  wing  in  Figure  1(c). 


(a)  (b)  (c) 

Figure  1.  Three  representative  MAVs.  (a)  Fixed  wing*21;  (b)  rotary  wing*31;  (c)  flapping  wing**l 

Most  of  the  published  research  on  wind  gust  is  on  the  dynamic  loading  and  structural  safety  of  manned 
aircraft,  which  means  that  the  Reynolds  number  is  in  much  higher  regimes.  For  example,  Zaide  and  Raveb 
used  a  Navier-Stokes  equation  solver  to  simulate  the  response  of  an  airfoil  to  arbitrary-shaped  gust  * 1,1  with 
Mach  number  from  0.1 1  to  0,7.  Patil  and  Taylor  calculated  the  flight  dynamic  characteristics  and  gust 
response  of  highly  flexible  aircraft.*1’1  Golubev  et  al.  conducted  parametric  study  of  nonlinear  gusl-airfoil 
interaction  for  a  range  of  gust  frequencies  and  amplitudes  in  subsonic  flow.*131  Yang  and  Obayashi  coupled 
a  Navier-Stokes  equations  and  dynamic  equations  to  study  the  gust  responses  in  plunging  and  pitching  for  a 
supersonic  transport  vehicle  with  and  without  the  consideration  of  structural  deformation.1'"*  On  the  MAV 
front,  Shyy  et  al.  computationally  and  experimentally  investigated  the  aerodynamic  response  of  rigid  and 
flexible  airfoils  to  gust*1*1  and  reported  that  a  flexible  airfoil  is  more  gust-tolerant  than  a  rigid  airfoil  in 
terms  of  maintaining  the  lift  to  drag  ratio.  In  their  study  of  a  low  Reynolds  number  airfoil  SD7003  in  the 
transitional  flow  region,  Ltan  and  Shyy1*'1*  found  that  both  the  lift  and  drag  coefficients  have  substantial 
variations  under  gust  environment.  The  gust  causes  the  formation  and  burst  of  the  laminar  separation 
hubble  and  ultimately  affects  the  airfoil  performance. 
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To  support  the  design  of  MA  Vs,  we  endeavor  to  understand  the  unsteady  aerodynamics  inherent  in  the 
small  flight  vehicles.  In  this  paper  we  use  computational  fluid  dynamics  (CFD)  to  study  the  fluid  physics 
of  Happing  airfoils  in  flow  regimes  relevant  to  MAVs.  Specifically,  our  goals  are  to 

(i)  qualify  the  role  played  by  Sirouhal  number,  and  plunging  and  pitching  frequency  and 
amplitude  in  the  How  structure  and  the  force  generation; 

(ii)  investigate  the  response  of  plunging  airfoil  under  both  steady  and  unsteady  freestream 
conditions; 

In  the  following,  we  first  present  some  background  in  the  aerodynamics  of  flapping  airfoil.  Then  we 
introduce  the  numerical  tools  used  in  the  work.  In  the  numerical  result  part,  we  validate  our  analysis  tool 
against  some  well  know  lest  cases.  At  last  we  compare  the  performance  of  flapping  wing  and  fixed  wing 
under  gust  environment. 


Flapping  Airfoil  Aerodynamics 

Flapping  Kinematics 

The  kinematics  of  an  airfoil  experiencing  a  combination  of  harmonic  plunging  and  pitching  motion  can  be 
describes  as  follows: 

h(t)/c  =  h0  cos (d)i  +  )  0(t)  -  #0  coster  +  )  ( 1 ) 

where  h(i  is  the  plunging  amplitude  non-dimensionalized  by  the  airfoil  chord  c,  depicts  the  pitching 
motion  of  an  amplitude  of  On  and  a  frequency  of  to.  For  the  flapping  motion  hn  is  defined  positive  upward 
and  pitch  amplitude  positive  clockwise.  The  leading  angle  between  pitch  and  plunge  is  defined  as 
V*  =  $0  ~$h .  The  kinematics  can  be  schematically  illustrated  in  Figure  2 


Figure  2.  Schematic  of  Happing  motion. 


In  the  Happing  wing  study,  two  non-dimensional  parameters  have  been  employed  in  the  literature:  one  is 
i he  reduced  frequency  defined  as: 


k  -  (0C 

another  is  the  Strouhal  number  based  on  Anderson  et  aL|7S: 

2  fcK 


Si 


(2) 


(3) 


where  Ult  is  the  reference  velocity,  typically  the  freestream  velocity  or  forward  flight  velocity.  Since  there 
are  two  length  scales,  chord,  cf  and  stroke  amplitude  ch0.  Si  and  k  characterize  the  relative  time  and  length 
scales. 


Different  angles  of  attack 
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To  better  understand  the  flapping  wing  aerodynamics,  it  will  be  useful  to  look  at  several  different  angles  of 
attack  attributed  from  different  causes.  The  kinematic  angle  of  attack  is  relative  to  the  pivot  point  and  is 
defined  as  the  following: 

a(t)-  @{t)  +  fl(i)  =  &(t)+  ian"'(-  —  cos(a)t  +  $,)+  tan"1  (^5/ si n(£t>/  +  $k ))  (4) 

UH  di 

The  kinematic  angle  of  attack  consists  of  two  terms.  The  first  term  on  the  right  side  of  the  equation,  0  ,  is 
due  to  the  pitching  motion,  and  the  second  term,  0  ,  is  due  to  the  plunging  motion.  By  definition  the 
kinematic  angle  of  attack  is  exclusively  determined  by  the  freestream  velocity  and  the  airfoil  motion 
without  considering  the  wake  effect.  More  importantly,  it  neglects  the  time  history  of  the  flow.  As  we  will 
discover  later,  flow  history  leads  to  phase  lag  in  the  force.  From  Eq.  (4)  we  can  readily  see  that  the 
Strouhal  number  affects  the  kinematic  angle  of  attack.  Its  role  is  demonstrated  in  Figure  3.  At  a  small 
Strouhal  number,  the  contribution  from  the  plunging  motion  is  small,  and  the  kinematic  angle  of  attack 
closely  follows  the  shape  of  a  sinusoidal  function.  As  the  Strouhal  number  increases,  the  kinematic  angle 
of  attack  no  longer  follows  the  sinusoidal  shape. 


Figure  3.  An  instantaneous  kinematic  angle  of  attack  at  different  Strouhal  numbers.  The  angle  of  attack 
becomes  less  sinusoidal  as  the  Strouhal  number  increases.  (h0=  0.75,  and  the  nominal  angle  of  attack, 
a0  =-0p  +  tan ‘‘(2*^),  is  15°) 

Thrust  by  lift 

It  is  known  that  a  flapping  airfoil  can  generate  thrust  at  certain  combinations  of  oscillation  frequency  and 
amplitude.  This  problem  has  been  extensively  investigated  experimentally  ,7J7  I*J9^  and  numerically 
typically  by  visualizing  the  wake  structures.  Without  invoking  sophisticated  analysis  tools,  we  can  get  an 
intuitive  idea  of  how  the  thrust  is  generated  by  looking  at  the  local  flow  direction,  a  concept  frequently  used 
in  the  aerodynamics  textbooks  to  calculate  the  induced  drag.  Consider  a  flapping  airfoil  as  shown  in  Figure 
4.  The  naked  eye  sees  the  airfoil  at  an  angle  of  attack  of  aK  ,  called  geometric  angle  of  attack  and  defined 

between  the  mean  chord  and  the  freestream  flow  direction.1211  However,  the  Happing  motion  induces  a  flow 
velocity  Uf  near  the  wall.  Due  to  the  flapping  motion,  the  airfoil  actually  encounters  the  deflected  local 
flow  Uf,  which  is  the  superposition  of  the  freestream  velocity  U  and  the  induced  velocity  Ur  During  the 
plunging  down  motion,  flow  in  the  vicinity  of  airfoil  is  inclined  upward  with  respect  to  the  freestream,  and 
the  lift  vector  that  remains  perpendicular  to  the  local  flow  direction  is  liked  forward.  And  this  tilted-lift 
contributes  to  thrust.  Similarly,  flow  is  inclined  downward  during  the  plunging  up  motion  and  the  lift 
(negative)  induces  thrust  too,  Theodorsen1221  and  Garrick,1*31  provided  analytical  expressions  for  lime 
dependent  lift,  moment,  and  thrust  on  a  sinusoidally  plunging  and  pitching  flat  plate  in  a  2D  potential  flow 
Their  analysis,  though  has  limitations,  is  a  good  starting  point  to  study  the  flapping  airfoil.  It  should  be 
noted  that  Figure  4  and  the  associated  analysis  facilitate  our  understanding  of  thrust  generation  but  to 


4 


calculate  the  thrust  requires  sophisticated  analysis  tools  with  the  consideration  of  viscous  force  and 
unsteady  aerodynamics. 


f 


Induced  wi 


Chord  line  direction 


Lorn  I  lluiv  dherUon 


Figure  4,  The  origin  of  thrust  from  a  flapping  airfoil. 


During  the  flapping  motion,  the  airfoil  experiences  the  instant  force  per  unit  span  in  the  forward  direction, 
force  in  the  transverse  direction,  pitching  moment  around  pivot  point,  and  instantaneous  power  input, 
which  arc  notated  as  7(7),  L(t)r  A/(t),  and  P(i),  respectively.  Among  those,  the  input  power  is  computed  as 
follows: 


dt  dt 

The  thrust  coefficient  Cj  and  the  power  input  coefficient  Cp  are  normalized  as  followed: 


Ct  =  i 


TO) 


-pvic 


c-i 


p(0 


-pule 


2  2 
The  lime-averaged  thrust  coefficient  CTmfm  and  input  power  coefficient  Cpmtm  can  be  evaluated  as: 


^T  mrtyt  “ 


=  rf  C,  (»di 


where  T  =  Ini  <a  .  Finally,  the  propulsive  efficiency,  >/,  is  defined  as: 

r 

7  mum 

n^- - 


(5) 

(6) 

(7) 

(8) 


Reduced  frequency  and  Slrouhal  number 

There  have  been  different  opinions  on  the  roles  of  these  two  parameters.  Taylor  et  al.  examined  existing 
data  of  wingbeat  frequency,  stroke  amplitude  and  flight  speed  of  42  species  from  bird,  bat,  and  insect 
families.  They  found  that  75%  of  the  species  falling  in  the  range  of  0.19<5/<0.41  for  cruising  flight,  with 
an  optimum  of  around  0.3. [241  Given  the  disparate  variety  of  species  used  in  the  study,  they  implied  the 
Slrouhal  number  was  turned  for  high  propulsive  efficiency  for  cruising  flight  or  swimming.  Experiments 
with  isolated  airfoils  also  reported  high  propulsive  efficiencies  within  the  interval  of  [0.2  0.4] 
Triamafyllou  ct  al.  stated  that  the  Slrouhal  number  almost  solely  determined  the  thrust  efficiency.^51 
Experiments  by  Anderson  et  al.  [7Kvith  a  NACA0012  airfoil  undergoing  combined  plunging  and  pitching 
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motions  at  Re=4*I04  also  showed  that  the  efficiency  peaked  at  the  Strouhal  number  range  0.25<S><0,4, 
with  an  high  efficiency  of  87%/3  These  evidence  support  the  statement  that  the  Strouhal  number  is  the 
dominant  parameter  in  using  flapping  wing  for  efficient  propulsion,1  201 

On  the  other  hand,  Bandyopadhyay  el  al.  performed  experiments  by  attaching  a  pair  of  flapping  foils  to  the 
rear  of  a  rigid  cylinder  to  mimic  the  fish  swimming.*271  In  their  experiment,  alternate  vortex  was  shed  from 
the  nose  of  the  cylinder.  They  found  that  the  Strouhal  number  was  by  no  means  the  only  one  to  determine 
the  thrust  efficiency  when  there  was  an  interaction  between  vortices  shed  from  the  leading-edge  and 
flapping  foils.  They  observed  that  the  peak  efficiency  was  reached  below  the  range  0.25<5if<035.  In 
another  work,  Triantafyllou  el  al.  suggested  that  both  the  flapping  frequency  (hence  the  reduce  frequency) 
and  the  Strouhal  number  were  needed  to  describe  the  flow.12*1  Experimental  study  by  Lai  and  Platzer,y|  and 
numerical  simulation  by  Young*2”1  show  that  both  the  reduced  plunging  frequency  k  and  non- 
dimensional  i  zed  amplitude  /j0  should  be  considered  separately  to  determine  the  wake  structure. 

From  our  early  discussion  we  know  that  the  reduced  frequency  characterizes  the  temporal  property  of  the 
flapping  motion  while  the  Strouhal  number  considers  the  combined  effect  from  temporal  and  spatial 
perspectives  but  ignores  their  individual  effect.  In  the  authors1  belief,  both  temporal  and  spatial  properties 
should  be  individually  considered  in  the  flapping  wing  study.  This  can  be  understood  from  the  following 
hypothetical  example.  Suppose  the  Strouhal  number  is  the  sole  parameter  determining  the  propulsive 
efficiency,  from  Eq.  (3)  we  know  that  for  airfoil  with  chord  c  and  forward  speed  of  UUy  as  long  as  they?*,,  is 
the  same,  it  has  the  same  efficiency.  If  one  combination  of  frequency /and  hu  produces  considerably  high 
efficiency  then  we  can  always  increase  the  value  of  h0  and  correspondingly  decrease  /  to  make  jhn 
constant.  Under  extreme  situation  we  can  make  h0  infinitely  large  and/ infinitesimally  small  so  that  the 
airfoil  is  essentially  stationary.  Because  a  stationary  airfoil  does  not  generate  thrust,  we  know  the  argument 
that  the  Strouhal  number  is  the  sole  parameter  determining  the  efficiency  is  not  true. 

Numerical  Methods 


Reynolds- Averaged  Navicr-Stokes  Solver 

The  full  Navier-Stokes  equations  in  curvilinear  coordinates  is  solved  with  a  pres  sure- based  algorithm, 
generalizing  from  the  original  Semi- Implicit  Method  for  Pressure- Linked  Equations  (S1MPLE).*wjoJ  The 
convection  terms  are  discretized  with  the  second-order  upwind  scheme  and  the  diffusion  terms  are 
discretized  with  the  second-order  central  difference  scheme,  The  time  integration  is  performed  with  a 
second-order  implicit  three-point  backward  scheme  for  better  handling  of  accuracy  and  strict  time  step 
constraint  imposed  by  the  extremely  fine  grid  resolution.  The  turbulence  computation  is  realized  with 
Wilcox’s  k-to  turbulence  mode  I1 31  ]  and  the  transitional  How  simulation  is  carried  out  through  a  method  by 
Lian  and  Shyy,***1  which  augmented  a  low  Reynolds  number  k~m  turbulence  model*311  with  the  t?A  transition 
method.*321 

Particle  tracking  and  StreakJine  integration 

Lai  and  Platzcr  visualized  the  wake  structure  of  airfoil  N  AC  A  00 12  exerting  plunging  motion  by  releasing 
non- intrusive  dye  to  the  flow  from  the  injection  point.11”1  The  snapshot  of  the  wake  structure  is  the  path 
traveled  by  all  particles  that  have  passed  through  the  injection  point.  It  is  the  streakline  if  we  ignore  the 
mass  diffusion  of  the  dye.  In  the  numerical  analysis,  the  streakline  can  be  evaluated  with  the  following 
formula: 


*(*„,/)=  x„+ J  U(x,/)|*  (9) 

where  x0  is  the  injection  position,  u(x,t)  is  the  velocity  vector  at  x,  Equation  (9)  is  integrated  with  the  forth- 
order  Runge-Kutta  method.  The  velocity  vector  at  arbitrary  position  x  is  obtained  with  bilinear 
interpolation. 


Results  and  Discussion 
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For  demonstration  purpose,  we  use  the  NACA0012  airfoil  throughout  this  work.  The  readers  should  be 
aware  that  in  real  design,  a  slender  airfoil  with  modest  camber  is  more  suitable  for  MA  Vs  than  thicker  and 
symmetric  airfoils  like  NACA0012.  We  investigate  the  Happing  wing  aerodynamics  at  different  Reynolds 
numbers,  Strouhal  numbers,  reduced  frequencies,  and  freestream  environments.  The  Reynolds  number 
considered  in  this  paper  is  L2*10\  2.0*  10J  and  4.0* IQ4,  respectively.  The  Strouhal  number  covers  the 
range  from  0.03  and  0.5,  and  the  reduced  frequency  has  the  span  from  0,2  to  3.93.  Because  the  gust  in 
urban  environment  has  a  low  frequency  around  1  Hz,  the  study  of  flapping  airfoil  under  gust  condition 
involves  multiple  time  scales.  In  this  section,  we  first  validate  our  analysis  tools  and  then  investigate  the 
gust  impact  on  the  aerodynamics.  For  code  validation,  we  will  compare  the  vortex  shedding  frequency,  lift 
and  drag  history,  and  wake  structure  at  a  variety  of  plunging  amplitudes  and  frequencies.  Comparisons  are 
also  made  with  available  data  from  literature. 

Frequency  validation 

Koochesfahani  measured  the  natural  shedding  frequency  over  a  stationary  NACA0012  airfoil  at  Re  = 
1 .2*  1G'1.1  J  Even  though  the  airfoil  is  streamlized,  at  such  a  low  Reynolds  number,  it  acts  as  a  bluff  body 
due  to  the  separation  at  the  trailing  edge.  Vortices  are  shed  alternatively  from  both  sides  of  the  (railing 
edge,  Koochesfahani  reported  the  non-dimensionalized  natural  shedding  frequency  of  8,7.  With  a 
compressible  flow  solver  and  laminar  flow  simulation.  Young  obtained  a  frequency  of  8.5  at  a  Mach 
number  of  0.05  and  8.2  at  a  Mach  number  of  0.2, At  this  Reynolds  number,  flow  separates  at  the  trailing 
edge  without  undergoing  transition.  Our  simulation  shows  a  frequency  of  8.75.  Lai  and  Platzer1'*' 
visualized  the  wake  structure  over  the  airfoil  at  the  Reynolds  number  of  2*  JO4,  but  they  did  not  measure  the 
shedding  frequency.  Based  on  his  simulation  Young  reported  a  frequency  of  9.4.l:,}'  Our  laminar  flow 
simulation  reveals  that  the  natural  shedding  frequency  is  9,8.  In  neither  case,  flow  experiences  transition, 
hence  our  simulation  is  based  on  laminar  flow  assumption. 

The  approach  of  Young  followed  the  method  of  Tuncer  and  Platzer,'33'  in  which  the  Crank-Nicolson 
second-order  scheme  was  used  for  time  discretization  and  the  viscous  flux  terms  were  evaluated  using  the 
second-order  central  difference  scheme  and  the  inviscid  fluxes  with  a  third-order  accurate  Osher  upwind 
scheme.  However,  Tuncer  and  Platzer  used  the  thin-layer  approximation  while  Young  solved  the  full 
Navier-Stokes  equations.  Young  used  C-grid  with  541*61  points,  among  them  377  points  were  around  the 
airfoil  s ur face Jlr)f  The  first  normal  gridpoint  was  9,2*10'*  chord  away  from  ihe  surface.  In  our  computation, 
we  use  O-grid  with  401*250  points.  The  first  normal  gridpoinl  is  1.0*10*  chord  away  from  the  surface. 
Young  investigated  the  effect  of  freestream  Mach  number  on  the  numerical  results.^"  He  found  that 
increasing  the  frequency  or  increasing  the  Mach  number  decreased  the  wavelength.  As  the  wavelength 
decreased,  the  interaction  between  the  pressure  filed  immediately  close  to  the  airfoil  and  the  sequence  of 
pressure  wave  caused  by  the  previous  airfoil  motion  had  increasing  effect  on  changing  the  force  magnitude 
and  varying  the  phase  angles.'7"' 


Force  validation 

Lai  and  Platzer  experimentally  studied  the  wake  structure  based  on  a  NACA0012  airfoil  at  the  chord 
Reynolds  number  of  2*I0J  and  various  oscillation  frequencies  and  amplitudes.1'^  Corresponding  numerical 
simulations  were  performed  by  Young  and  Lai.1*1  To  validate  our  code  for  unsteady  force  computation,  we 
choose  two  tests.  In  the  first  case  the  reduced  frequency,  A,  is  3,93,  and  the  non-dimensional  plunging 
amplitude,  A0,  is  0.0125,  which  give  a  low  Strouhal  number  of  0.03,  Figure  5  shows  the  time  histories  of 
lift  and  drag  coefficients.  Same  as  flow  over  a  circular  cylinder,  the  drag  varies  with  twice  the  frequency  as 
lift-  At  this  frequency  and  amplitude,  the  airfoil  is  drag-producing.  Since  there  is  no  experimental 
measurement  for  force,  we  compare  our  results  with  the  numerical  results  by  Young.1341  Both  lift  and  drag 
coefficients  show  good  agreement  except  the  discrepancy  in  the  troughs  of  the  drag  coefficient  (Figure  5). 
With  the  implicit  scheme,  we  have  100  time  steps  per  flapping  cycle  while  Young  used  explicit  scheme 
with  16,860  time  steps  per  flapping  cycle. 
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Figure  5.  Drag  and  lift  coefficient  history  at  k  =  3.93,  h0  -  0.0125,  and  St  =  0.03. 

Force  probe  reveals  that  both  pressure  and  viscous  forces  contribute  to  the  total  drag  while  the  viscous 
force  contributes  almost  70%  of  the  total  drag.  The  time  averaged  viscous  drag  in  our  simulation  is  about 
4%  higher  than  that  of  Young,  At  this  Reynolds  number,  we  may  not  have  adequate  resolution  to  resolve 
the  flow  structure.  With  the  increase  of  plunging  amplitude,  the  portion  of  viscous  force  in  the  total  drag  is 
reduced.  As  shown  in  Figure  6,  the  difference  between  ours  and  that  by  Young  actually  decreases  at  higher 
oscillation  amplitude.  At  this  amplitude,  the  lime  average  pressure  force  is  thrust -generated. 


Q 

Q 


Figure  6.  Drag  and  lift  coefficient  history'  at  k  -  3.93,  ho  -  0.025,  and  St  -  0.06, 

As  discussed  in  Eq.  (4),  the  airfoil  experiences  an  instantaneous  kinematic  angle  of  attack  due  to  its 
flapping  motion.  Intuitively  the  higher  such  an  angle  of  attack  is,  the  higher  the  lift  coefficient  will  be. 
However,  this  intuition  is  based  on  two  assumptions:  quasi-steady  state  and  no  impact  from  the  wake. 
Figure  7  dearly  demonstrates  that  the  lift  variation  does  not  coincide  with  the  kinematic  angle  of  attack, 
with  the  lift  leading  the  angle  of  attack  by  around  40°.  Instead,  the  lift  closely  follows  the  airfoil  plunge 
acceleration,  d2h/dt 1 . 
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Since  the  pressure  force  dominates  the  total  lift,  it  will  be  valuable  to  examine  the  pressure  variation  during 
the  airfoil  plunging  motion.  Five  representative  time  instants  are  chosen,  as  marked  A,  B,  C,  D,  and  E  in 
Figure  7.  The  corresponding  pressure  coefficient  distributions  are  shown  in  Figure  8.  At  instant  A,  the 
airfoil  is  in  the  middle  plunging  plane  (hf(h{f:)  =  0)  and  both  the  lift  and  drag  coefficients  approaches  to 
their  minimum  magnitudes.  There  is  a  cross-over  point  in  the  pressure  coefficient  locating  at  around  35% 
of  the  chord  After  the  cross-over  point,  the  pressure  on  the  upper  surface  is  actually  larger  than  the 
pressure  at  the  corresponding  lower  surface,  which  reduces  the  total  lift.  During  the  plunging  down  motion 
(instants  B  and  C  in  Figure  8(a)),  the  cross-over  point  moves  toward  the  leading  edge,  which  generates 
negative  lift.  Once  the  airfoil  starts  to  plunge  upward,  the  cross-over  point  moves  toward  the  trailing  edge 
(instants  D  and  E  in  Figure  8(b)).  For  the  studied  case,  flow  history  is  more  important  than  the  kinematic 
angle  of  attack  in  terms  of  aerodynamic  forces,  especially  the  lift.  We  can  draw  the  same  conclusion  for 
other  tests  with  increased  plunging  amplitude  (0.05c  and  0.075c)  but  same  other  parameters.  These  tests 
suggest  that  quasi -steady  simulation  may  not  be  adequate  for  unsteady  simulation,  especially  when  the 
reduced  frequency  is  considerably  high.  As  we  will  see  soon,  in  the  combined  plunging  and  pitching  study 
with  plunging  amplitude  of  0.75c,  the  kinematic  angle  of  attack  seems  more  responsible  for  the  lift. 


Figure  7,  Drag  and  lift  coefficient  history  at  k  ~  3.93,  h()  =  O.OI25,  and  ,57  -  0.03. 

We  further  inspect  the  How  structure  relative  to  the  airfoil.  This  is  obtained  by  transforming  the  velocity 
field  to  the  coordinate  fixed  on  the  plunging  airfoil.  Figure  9  vividly  shows  the  pattern  of  the  kinematic 
angle  of  attack  changes  from  positive  to  negative,  as  is  consistent  with  figure  7.  Through  the  plunge  cycle, 
flow  separates  near  the  trailing  edge.  The  separated  flow  creates  a  blunt-edged  body  and  vortices  are  shed 
alternatively  from  the  trailing  edge.  A  close-up  reveals  that  there  are  two  vortices  near  the  trailing  edge  at 
instant  A  when  the  airfoil  is  at  its  middle  plane  and  starts  to  plunge  down.  As  shown  in  Figure  10  (a),  one 
vortex  rotates  clockwise  which  is  due  to  the  flow  separation,  and  another  vortex  is  small  and  rotates  anti¬ 
clockwise,  indicating  it  is  caused  by  the  airfoil  motion.  When  the  airfoil  passes  the  middle  plane  and 
moves  up,  the  small  vortex  rotates  clockwise  (Figure  JO  (b)). 

No  leading  edge  vortex  is  observed  when  the  plunging  amplitude  is  small  (0.0125  and  0.025).  But,  it 
appears  with  the  increase  of  plunging  amplitude.  Figure  1 1  shows  the  streamlines  at  the  relative  coordinate. 
Leading  edge  vortex  and  trailing  edge  vortex  coexist.  In  the  work  of  Young  no  leading  edge  vortex  was 
found  at  even  high  amplitude  of  0.25,  We  gather  that  the  leading  edge  vortex  is  caused  by  the  kinematic 
angle  of  attack  defined  in  Eq.  (4).  Hence  both  large  plunge  amplitude  and  frequency  will  contribute  to  the 
appearance  of  leading  edge  vortex  in  similar  way. 
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(a)  (b) 

Figure  8.  Pressure  coefficient  evolvement  during  half  plunging  cycle;  k  =  3,93,  /r0  -  0,0125,  and  57  =  0.03. 
(a)  One  quarter  cycle  from  instant  A  to  C  (from  middle  plane  to  bottom  plane);  (b)  another  quarter  cycle 
from  bottom  plane  to  middle  plane. 
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(b) 

Figure  10,  Close-up  of  instantaneous  streamlines  in  the  relative  coordinate  during  the  half  flapping  cycle. 
At  k  =  3.93,  h0  =  0.0125,  and  Si  -  0,03.  (a)  Plunge  down  at  hl(hac)= 0;  (b)  Plunge  up  at  A/(Aac)=0. 


Figure  1 1.  Instantaneous  streamlines  in  the  relative  coordinate  when  the  airfoil  is  at  the  middle  plane  and 
pi  unges  down .  k  -  3 .93,  h0  -  0 .05,  and  St  =  0 . 1 25 . 

Jones  et  al.  experimentally  visualized  the  wake  structures  of  a  plunging  airfoil.  Depending  on  the  flapping 
amplitude  and  frequency,  they  found  the  drag-producing  Karman  street,  the  neutral  wake  and  the  thrust- 
producing  reverse  Karman  street.  They  indicated  that  the  wake  structure  depended  on  the  kh0  (=2nSt), 
However,  experiments  by  Triantafyliou  et  Ohashi  and  IshikawaJ3^  and  Kadlec  and  Davis1**1  showed 
the  wake  structure  was  independently  determined  by  the  both  k  and  h 0>  which  was  further  confirmed 
experimentally1  iyJ  and  numerically.1^1 
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Figure  12  illustrates  the  streak! ines  in  the  wake  region  at  two  different  flapping  amplitudes*  both  of  which 
show  the  reverse  Kantian  street  and  indicate  thrust  generation.  For  the  considered  two  cases,  flow 
visualization  reveals  a  clear  wake  structure  near  the  trailing  edge  j 171  especially  indicating  that  the 
diffusivity  rate  at  that  region  is  low.  We  assume  flow  is  laminar.  Simulations  by  Young  showed  little 
difference  between  laminar  and  turbulent  results,1  :"3  A  large  scale  similarity  with  the  experimental 
visualization  is  apparent. 


CFD:  Lian  and  Shyy: 
A  “  3,93*  h0=  0.05, 
and  St  =  0.125, 


CFD:  Lian  and  Shyy: 
A  =  3.93,  A*-  0,075, 
and  St  =  0.1 88. 


Figure  12.  Streak  line  plot  of  flow  over  plunging  airfoil  NACA0012  at  Re  =  20,000  and  reduced  frequency 


of  3.93. 


Combined  plunging  and  pitching  motion 

Anderson  et  al.  measured  the  time-averaged  thrust  coefficient,  input  power  coefficient,  and  propulsion 
efficiency  of  a  NACA0012  airfoil  undergoing  combined  sinusoidal  plunging  and  pitching  motion  in  the 
testing  tank  facility  at  The  tank  has  the  dimensions  of  30  m  *  2.5  m  *  1.2  m  and  can  run  at  constant 

speed  up  to  3  m/s.  The  airfoil  has  a  chord  of  10  cm  and  span  of  60  cm.  Circular  end  plates  are  fined  to 
avoid  the  three-dimensional  end  effects,  A  piezoelectric  force  transducer  provides  force  measurement  in 
three  axes  while  another  transducer  measures  the  force  transmitted  to  the  chain  to  find  torque.  The  reported 
Reynolds  number  is  Re  =  4*  104,  and  the  pitch  axis  is  located  at  the  1/3  chord  point. 
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Among  their  systematic  tests  with  varied  plunging  amplitude,  nominal  angle  of  attack,  and  Strouhal 
number,  one  of  the  optimal  eases  is  obtained  for  the  plunging  amplitude-to-chord  ratio,  =  075,  the  phase 
angle  between  plunging  and  pitch  (pitch  leading  plunge)  is  75°,  and  the  normal  angle  of  attack,  a®  =  15°. 
We  study  is  centered  with  this  setting  with  varied  Strouhal  number.  The  readers  should  be  aware  that  based 
on  full  turbulent  simulation  with  the  k-to  turbulence  model,  the  turbulent  Reynolds  number  in  the  vicinity 
of  the  airfoil  is  in  the  order  of  100,  indicating  the  Reynolds  number  of  4*  IG4  may  not  be  high  enough  to 
prompt  full  turbulent.  On  the  other  hand,  laminar  flow  simulation  may  not  have  the  resolution  to  resolve 
the  flow  structure  at  this  Reynolds  number.  We  do  not  claim  that  we  fully  resolve  the  low  Reynolds 
number  aerodynamics.  Rut,  our  main  purpose  is  to  look  at  the  aerodynamics  based  on  available  analysis 
tool  For  the  flapping  wing  study,  we  use  both  laminar  and  turbulent  simulations  but  the  analysis  except 
that  shown  in  Figure  13  is  based  on  laminar  flow  simulation. 


(a)  (b) 

Figure  13.  (a)  Thrust  coefficient;  (b)  power  efficient  as  a  function  of  Strouhal  number  for  ho=  075,  = 

1 5°,  and  ^  =  75°. 

We  compare  the  mean  thrust  coefficient  and  input  power  at  different  Strouhal  number  with  both 
experimental  and  numerical  data  in  Figure  13.  For  all  the  Strouhal  number  except  at  Strouhal  number  of 
0.1,  the  laminar  and  turbulent  simulations  under-predict  the  mean  thrust  coefficient.  However,  the  input 
power  coefficients  from  our  laminar  flow  simulation  closely  match  the  measurement  by  Anderson  et  alJ7* 
Overall,  all  the  data  shows  that  the  lime-averaged  thrust  and  input  power  increase  with  the  Strouhal 
number. 

Figure  14  shows  the  time  histories  of  lift,  thrust,  input  power,  and  moment  during  one  flapping  cycle, 
starting  from  the  middle  stroke  plane.  The  Strouhal  number  is  03  and  reduced  frequency  is  0.63.  At  the 
beginning  of  the  cycle,  the  airfoil  plunges  down  (Figure  15(a))  with  the  maximum  negative  pitching  angle 
6  (Figure  15(b)),  However,  both  the  lift  and  thrust  are  near  their  maximum  values.  The  reason,  as 
explained  before,  is  that  the  airfoil  actually  experiences  a  large  kinematic  angle  of  attack,  a  (Figure  15(c)). 
From  Figure  15(d)  we  know  that  at  instant  A  the  angle  of  attack  due  to  plunging  motion,  p,  reaches  its 
maximum.  As  illustrated  in  Figure  16  (a),  the  local  flow  is  deflected.  The  focal  flow  strikes  the  airfoil 
from  below  instead  of  in  front.  The  lift,  which  is  perpendicular  to  the  local  flow,  is  then  inclined  toward 
upstream,  which  leads  to  thrust  (Figure  14(b)).  Similar  conclusions  can  be  drawn  at  instant  B. 

When  the  airfoil  reaches  its  lowest  position  during  the  cycle  (point  C),  the  kinematic  angle  of  attack 
becomes  negative.  From  Figure  16(c)  we  also  see  that  the  local  flow  meets  the  airfoil  from  above. 
However,  both  the  lift  and  thrust  coefficients  are  close  to  zero.  As  we  mentioned  before,  the  flow  history 
may  be  responsible  for  that.  After  point  C,  the  airfoil  starts  to  reverse  its  plunging  direction.  Meanwhile, 
the  leading  edge  of  the  airfoil  pitches  up  and  the  trailing  edge  pitches  down.  At  instant  D,  the  kinematic 
angle  of  attack  is  negative,  and  How  is  deflected  downward.  Meanwhile,  the  lift  is  tilted  towards  upstream. 
Again,  the  tilted  lift  leads  to  thrust.  It  is  a  half-cycle  from  point  A  to  point  D,  and  the  other  half-cycle  from 
point  E  to  point  H  can  be  seen  as  the  opposite  side  of  the  first  half  cycle. 
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Figure  14.  Time  histories  of  lift,  thrust,  moment,  and  input  power  coefficient  for  at  h0  =  0.75,  St  =  0J,  k  . 
0.63,  t//=  75°,  and  =  15°  using  laminar  simulation. 

It  is  instructive  to  fook  at  the  flow  filed  to  understand  how  the  leading  edge  vortex  enhances  the  lift  and 
thrust.  During  the  plunging  down  motion,  a  leading  edge  vortex  is  formed  and  travels  downstream  (Figure 
16),  The  vortex  core  is  usually  associated  with  low  pressure  (Figure  17),  which  enhances  the  lift  and  thrust. 
When  it  travels  downstream,  the  vortex  gradually  loses  its  coherent  structure  and  strength,  and  the  pressure 
approaches  its  ambient  value.  For  example,  the  pressure  within  the  core  at  instant  C  is  about  one  third  of 
the  value  at  instant  A.  The  vorticity  contours  at  the  corresponding  time  instants  are  shown  in  Figure  18, 
The  vortex  leaves  off  the  trailing  edge  and  is  shed  into  the  wake. 
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(C)  (d) 


Figure  15,  Time  history  for  a(  ha-  0.75,  Si  -  0.3,  k  -  0.63,  tir=75°,  and  o^=  15°. 


(a>  (b) 


Figure  16.  Streamlines  seen  from  the  coordinates  fixed  on  the  airfoil,  (a)  h!(hgc)  =  0;  (b)  h/(hac)  =  -0.6;  (c) 
A/(M  =  -1;  (d)  h/(hoC)  =  -0.4.  h0  =  0.75,  Si  =  0.3,  k  =  0.63,  v=  75°,  and  »„=  15°. 
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Figure  17.  Pressure  coefficient  distribution  at  three  time  instants,  h0  =  0.75,  Sf  =  03,  k  =  0,63,  if/-  75°,  and 
<*,=  15°. 


(c)  (d) 

Figure  18,  Vortidty  contours  at  different  instants  during  half  flapping  cycle,  (a)  hf{h{!c)  =  0;  (b)  k{(hifc)  -  - 
0.6;  (c)  h/(hfJc)  -  -\ ;  |d)  hf{ht!c)  =  - 0A ,  A0  =  0.75,  Si  =  03,  k  =  0.63,  if/  -  75°,  and  15°. 

Stationary  airfoil  in  unsteady  freestream 

We  study  the  response  of  a  stationary  NACA0012  m  unsteady  flow.  The  Reynolds  number  is  40,000 
based  on  the  frees  tream  velocity  and  chord  length  and  the  angle  of  attack  is  set  to  4°.  At  this  Reynolds 
number  and  angle  of  attack,  laminar  separation  bubble  causes  flow  to  experience  laminar  to  turbulent 
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transition.  We  use  a  method  by  Li  an  and  Shyy^  to  simulate  the  transitional  flow.  Figure  19  shows  the 
computed  streamlines  and  the  normalized  turbulence  shear  stress  ,  Our  simulation  predicts  that 

transition  occurs  at  65%  of  the  chord  based  on  a  relatively  low  freestream  intensity  of  0.1%  while 
XFoil13*1  predicts  transition  position  is  66.7%.  The  lift  and  drag  coefficients  are  0.53  and  0.03  from  our 
simulation  and  0.53  and  0.026  from  XFoiL  respectively. 


We  limit  our  focus  on  head-on  gust  with  a  single  frequency  of  I  Hz: 

U(t)-  £/„(!  +  NA  stn(arf))  (10) 

it  varies  sinusoidally  with  time  with  a  frequency  of  m  and  variation  of  NA.  Here  m  is  equal  to  2rr  and  NA  is 
0.2.  Figure  20(a)  compare  both  the  lift  and  drag  coefficients.  The  comparison  between  the  unsteady  and 
quasi-steady  is  made  at  the  same  Rey  nolds  number.  Overall,  the  quasi -steady  simulation  predicts  larger  lift 
and  lower  drag  than  the  unsteady  simulation.  The  lift  curve  of  the  quasi-steady  simulation  has  similar 
pattern  as  the  unsteady  computation,  both  following  the  variation  of  U 2  /£7* .  On  the  other  hand,  ihc  drag 
curve  from  the  unsteady  simulation  is  quite  different  from  the  quasi-steady  prediction. 


(a)  (b) 

Figure  20.  Comparison  of  lift  and  drag  coefficients  between  unsteady  and  quasi-steady  computations  of  a 
stationary  airfoil  (a):  lift  coefficient;  (b):  drag  coefficient. 

Flapping  airfoil  in  unsteady  freestream 

The  flapping  airfoil  follows  plunging  motion  described  previously  in  Eq.  (1).  The  pitching  motion  is  as 
follows: 


0(f  1  =  Oti  cos(  cat  +  ^ )  +  4'1 


(ID 
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The  other  parameters  are  set  to  Si  ~  03,  k  -  0.63,  h0—  035,  =  1 5°,  co  -  2k,  and  y  -  75°.  In  the  case,  the 

time  averaged  pitching  angle  of  attack  is  4°.  The  same  head-on  gust  in  Eq.  (10)  is  applied  to  the 
freestream.  Hence,  the  gust  has  a  frequency  ofl  Hz  while  the  airfoil  flapping  at  7.4  Hz. 

Time  histories  of  the  lift  and  thrust  are  shown  in  Figure  21,  in  which  both  the  lift  and  thrust  are 
nondimensionalized  by  their  mean  values.  During  each  flapping  cycle,  neither  the  lift  nor  the  thrust  has 
the  symmetric  pattern  as  shown  early  in  Figure  14.  We  average  the  lift  and  thrust  in  each  flapping  cycle 
and  plot  them  in  Figure  22,  None  follows  the  variation  of  U1  fU\  .  Instead,  the  time-averaged  thrust 
decreases  as  the  freestream  velocity  increases.  It  seems  there  is  a  phase  lag  between  the  gust  variation 
and  the  response  from  the  airfoil. 


(a)  (b) 

Figure  21.  Force  history  of  a  flapping  airfoil  during  one  gust  cycle,  (a)  Lift  (b)  thrust.  h0  =  0.75,  Si  -  0.3,  k 
=  0.63,  y/  75°,  and  0^=  15°.  The  ratio  between  the  flapping  frequency  and  gust  frequency  is  7.4 


Figure  22.  Time-averaged  force  of  a  flapping  airfoil  during  one  gust  cy  cle,  (a)  Lift  (b)  thrust,  h0  =  0.75,  Si  = 
0.3,  k  =  0.63,  75°,  and  1 5°* 


Conclusion 


We  use  the  NACA0012  airfoil  as  one  example  to  study  the  flapping  airfoil  aerodynamics.  The  Reynolds 
number  is  in  the  range  between  10J  and  4.0*  104;  the  non-dimensional  Happing  amplitude  spans  from 
0.0125  to  0.75;  the  Strouhal  number  covers  a  range  between  0.03  and  0.5.  Based  on  our  theoretical 
analysis  and  numerical  simulation  we  have  the  following  findings: 
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])  Both  the  reduced  frequency  and  the  Strouhal  number,  characterizing  the  time  scale  between  flapping 
and  free st ream  motion  and  length  scales  between  chord  and  stroke  amplitude,  should  be  considered 
independently  to  determine  the  Happing  airfoil  aerodynamics,  which  includes  the  wake  structure,  lift 
and  thrust  generation,  and  leading  and  trailing  edge  vortex  formation. 

2)  The  plunging  airfoil  can  generate  thrust  or  drag  depending  on  the  plunging  frequency  and  amplitude. 
At  fixed  flapping  frequency,  increasing  the  plunging  amplitude  can  make  the  airfoil  from  drag 
generating  to  thrust  generating. 

3)  When  the  flapping  frequency  is  fixed,  at  small  the  flapping  amplitudes  ofG.0125  and  0.025),  there 
is  only  trailing  edge  vortex  due  to  the  effective  bluff  trailing  edge  and  the  plunging  motion;  at  large 
amplitudes  (hn  of  0.05  and  0.075),  there  are  both  leading  edge  vortex  and  trailing  edge  vortex. 

4)  When  the  reduced  frequency  is  considerably  high  (k  of  3.93),  the  time  history  effect  has  profound 
impact,  and  the  lift  variation  does  not  follow  the  change  of  the  kinematic  angle  of  attack. 

5)  When  the  plunging  amplitude  is  fixed,  increasing  the  Strouhal  number  from  0 .1  to  0.5  will  increase 
the  thrust  and  input  power  coefficients. 

6)  During  the  flapping,  a  leading  edge  vortex  enhances  the  thrust  because  it  generates  a  low  pressure 
zone.  This  vortex  gradually  loses  its  strength  when  it  travels  downstream. 

7)  The  response  of  a  stationary  airfoil  in  gust  environment  can  be  approximated  with  quasi-steady 
approach  when  the  gust  varies  with  low  frequency. 

8)  The  lift  and  thrust  of  a  flapping  airfoil  do  not  follow  the  change  of  the  gust. 
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